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ABSTRACT 

Molecular  dynamics  (MD)  implemented  on  parallel  processors  was  used  to  model 
supercritical  droplet  phenomena  occurring  in  combustion  devices.  The  use  of  molecular  dynamics 
allows  the  modeling  of  supercritical  phenomena  without  an  a  priori  knowledge  of  the  equation  of 
state  or  transport  properties  of  the  individual  components  or  the  mixture.  Three-dimensional 
supercritical  oxygen  vaporization  into  gaseous  oxygen  and  helium  using  two-site  Lennard- Jones 
potentials  for  the  oxygen  has  been  modeled  and  both  the  disappearance  of  surface  tension  above 
the  critical  point  and  the  modification  of  the  critical  point  for  a  binary  mixture  have  been  observed. 
A  distinct  change  in  droplet  morphology  was  observed  when  passing  through  its  critical  point. 
The  droplet  remains  spherical  as  it  vaporizes  under  subcritical  conditions  but  becomes  broken  and 
cloud-like  when  supercritical.  Equations  of  state  and  transport  coefficients  for  mass,  momentum 
and  energy  have  been  calculated  for  supercritical  argon,  nitrogen  and  oxygen  which  agree  with 
NIST  values. 


RESULTS 

Argon  and  LOX  Evaporation.  The  results  were  obtained  on  either  8  or  32  nodes  of  IBM’s 
Scaleable  Powerparallel  2  (SP-2)  using  the  Message  Passing  Interface  subroutine  library  for 
parallel  communication.  In  order  to  parallelize  the  problem,  the  molecules  were  distributed  evenly 
over  all  of  the  processors,  regardless  of  their  location  in  the  computational  domain.  This 
technique,  called  ‘atom  decomposition’  [1],  allowed  almost  perfect  load-balancing  to  be  achieved. 
The  results  presented  here  use  a  common  pair-wise  additive  intermolecular  potential  called  the 
Lennard-Jones  12-6  potential  (Fig.  1)  [2].  This  potential  contains  two  parameters  that  determine 
the  molecular  type:  a,  the  zero  energy  separation  distance,  and  e,  the  minimum  energy.  The 
gradient  of  this  potential  is  used  to  determine  the  force,  which  is  then  summed  over  alt  pairs  to 
determine  the  net  force  on  a  given  atom.  This  is  applied  to  each  atom  of  an  oxygen  molecule, 
which  is  called  the  Lennard-Jones  site-site  approximation. 

5,500,  27,000  and  100,500  (366,800  system  total)  atom  liquid  argon  droplets  have  been 
evaporated  under  both  subcritical  and  supercritical  conditions  [3-5].  Evaporation  is  carried  out  by 
placing  an  equilibrated  droplet  within  a  warm  gaseous  environment.  For  the  subcritical 
simulations,  the  d^  versus  time  law  was  obtained.  Rgure  2  shows  the  change  in  surface  area  with 
time  for  two  subcritical  conditions  and  one  supercritical  condition  with  two  initial  droplet  sizes.  As 
expected,  as  the  environment  pressure  and  temperature  increase,  so  does  the  droplet  evaporation 
rate.  For  the  supercritical  simulations,  both  the  5,500,  27,000  and  100,500  atom  drops  reached 
the  same  surface  regression  rate  after  an  initial  heat-up  period.  Figure  3  shows  the  supercritical 


surface  regression  rates  for  three  different  size  droplets  scaled  by  the  initial  surface  area.  The 
27,109  and  100,570  atom  droplets  show  identical  behavior  when  scaled,  with  the  5,587  atom  drop 
following  the  same  curve  but  with  slightly  greater  statistical  variation.  This  leads  us  to  believe  that 
extremely  large  systems  are  not  required  to  accurately  model  phase  change  phenomena.  Also 
significant  is  the  result  that  the  surface  tension  disappeared  after  the  same  time  interval  for  the 
differently  sized  supercritical  drops.  Subcritical  droplets  retain  a  spherical  geometry  while 
supercritical  droplets  quickly  lose  their  spherical  shape  and  become  more  cloud-like  in  form. 

Although  Lennard- Jones  is  sufficient  to  describe  intermolecular  interactions,  some  other 
method  must  account  for  constraining  two  oxygen  atoms  to  be  in  the  same  molecule.  Ideally,  one 
would  like  to  treat  chemical  bonds  as  additional  terms  in  the  potential  energy  equation,  but  this 
proves  to  be  both  intractable  and  unnecessary.  Because  bond  vibrations  tend  to  be  both  high 
frequency  and  low  amplitude,  one  can  simply  fix  the  interatomic  distance  to  a  prescribed  value  with 
little  consequence  to  the  results  of  the  simulation.  This  approach,  of  course,  would  not  be  valid  for 
large  molecules  where  torsional  motion  about  the  bonds  may  have  to  be  included.  Ano^er 
difficulty  encountered  when  simulating  multiple  chemical  species  is  determining  the  appropriate 
Lennard-Jones  parameters  for  interactions  between  unlike  molecules.  The  results  presented  here 
use  the  Lorentz-Berthelot  mixing  rules  [2].  Finite-differencing  of  the  equations  of  motion  was 
accomplished  using  a  modification  of  the  ‘velocity  Verlet’  algorithm  called  RATTLE  [6].  This 
algorithm  constrains  the  bond  length  of  the  oxygen  molecules  to  a  given  value  within  a  specified 
tolerance  through  an  iterative  process.  Both  the  positions  and  the  velocities  are  determined  at  the 

current  time  step  to  order  dt  ,  and  the  round-off  error  is  minimized. 

Two  representative  cases  of  an  oxygen  droplet  evaporating  into  a  helium  environment  are 
summarized  in  Table  1.  The  systems  consisted  of  roughly  20,000  moleculaes  and  the  initial  LOX 
droplet  diameter  was  7.7  nm.  As  a  reference,  the  critical  temperature  and  pressure  of  oxygen  is 
154.77  K  and  5.087  MPa,  respectively.  The  initial  condition  for  both  simulations  consisted  of 
placing  a  saturated  droplet  at  100  K  into  an  equilibrated  environment  at  the  various  conditions  listed 
m  the  table.  The  method  of  the  placing  the  droplet  into  the  environment  assured  that  the  initial 


Figure  1 .  Lennard-Jones  12-6  interatomic  potential  and  force. 


Figure  2.  Argon  droplet  surface  area  as  a  function  of  time  for  three  environmental 
conditions  and  two  drop  sizes  showing  linear  d^  versus  time  behavior. 
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Figure  3.  Supercritical  argon  droplet  scaled  surface  area  as  a  function  of  time  for  three 
initial  drop  sizes  showing  no  dependance  on  initial  droplet  size. 


surface  tension  of  the  droplet  was  retained  for  the  start  of  the  simulation.  The  simulation  geometry 
was  a  cube  with  periodic  boundaries  and  the  time  step  was  2.5  femtoseconds.  The  potential  cut¬ 
off  was  2.5  0,  which  refers  to  the  intermolecular  distance  at  which  the  potential  between  two 
atoms  is  so  small  as  to  be  neglected.  Energy  was  added  to  the  system  in  order  to  induce 
evaporation  by  ‘heating  the  boundaries’,  which  amounts  to  scaling  the  velocities  of  the  molecules 
in  the  boundary  region  of  the  cube  to  a  value  corresponding  to  a  prescribed  temperature. 

A  droplet  in  its  own  vapor  under  supercritical  conditions  shows  significant  differences  in 
evaporation  behavior,  such  as  immediately  vanishing  surface  tension  and  cloud-like  density 
profiles  [4].  The  results  of  extending  the  modeling  to  mixtures  is  shown  in  Figures  4  and  5. 
These  figures  contain  contour  plots  of  oxygen  mass  fraction,  density,  temperature,  and  the  average 
force  in  a  cubic  volume  for  two  different  helium  environments  at  times  of  250  psec  and  1 10  psec, 
respectively.  Because  surface  tension  is  the  macroscopic  manifestation  of  the  unequal  forces 
exerted  on  surface  molecules  in  a  liquid,  the  average  force  should  provide  an  approximation  of  the 
surface  tension  present  in  the  system.  Hence  a  large  surface  tension  would  be  indicated  by  a  bright 
ring.  The  environment  conditions  of  300  K  and  4  MPa  in  Figure  4  are  below  the  critical  pressure 
of  oxygen,  so  subcritical  behavior  is  expected.  This  is  observed  in  the  figure  by  a  spherical  droplet 
profile  and  strong  surface  tension.  In  contrast.  Figure  5  contains  results  when  the  environment 
pressure  and  temperature  are  at  300  K  and  20  MPa,  which  are  both  well  above  the  critical  point  of 
oxygen.  Although  one  might  expect  supercritical  behavior,  this  is  not  observed.  The  droplet 
profile  remains  spherical  and  surface  tension  is  retained  throughout  the  droplet  lifetime.  But 
experiments  have  shown  that  binary  mixtures  containing  helium  can  easily  have  critical  pressures 
over  100  MPa  [7]. 


Case 

Temperature 

Pressure 

1 

300  K 

4.0  MPa 

2 

300  K 

20.0  MPa 

Table  1:  Summary  of  thermodynamic  conditions  for  two  different  helium  environments  used  in  the 
LOX  droplet  evaporation  simulations.  The  critical  temperature  and  pressure  of  oxygen  is  154.77 
K  and  5.087  MPa,  respectively. 

Equation  of  state  and  transport  properties.  This  work  involved  the  determination  of  transport 
coefficients  and  pressures  of  supercritical  fluids  by  molecular  dynamics  (MD)  simulations  using 
the  Green-Kubo  formulae  [8]  and  the  virial  equation  of  state  [2],  respectively.  The  transport 
coefficients  that  are  of  interest  include  self-diffusion,  shear  viscosity,  and  thermal  conductivity. 
The  coefficients  are  computed  using  Green-Kubo  formulae  in  which  a  correlation  function  is 
integrated  over  time.  A  time  autocorrelation  function  measures  the  relation  between  the  values  of 
some  dynamic  quantity  at  two  different  times. 


C(0=  - J  A(ro)A(fo  +  t)dt^ 

(1) 

This  can  be  written  in  the  form 

(2) 

where  t©  and  t  denote  an  arbitrary  time  and  the  time  difference,  respectively,  and  the  brackets 
denote  an  autocorrelation. 

The  self  diffusion  process  is  the  mass  transport  associated  with  the  motion  of  a  single 
particle  in  a  fluid,  and  hence  characterizes  a  one-particle  property  of  the  system.  Self  diffusion 
may  be  described  by  the  mass  and  velocity  of  a  marked  particle  of  the  fluid.  For  the  self-diffusion, 
the  Green-Kubo  formula  contains  the  velocity  autocorrelation  function. 
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Figure  4:  Contour  plots  of  oxygen  mass  fraction,  density,  temperature,  and  the  average  force  after  250 
picoseconds  of  simulation  time.  The  helium  environment  is  at  300  K  and  4.0  MPa,  and  the  droplet  is 
exhibiting  subcritical  evaporation  behavior. 
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Figure  5:  Contour  plots  of  oxygen  mass  fraction,  density,  temperature,  and  the  average  force  after  110 
picoseconds  of  simulation  time.  The  helium  environment  is  at  300  K  and  20.0  MPa,  which  is  well  above  the 
critical  point  of  both  helium  and  oxygen.  The  droplet  is  still  exhibiting  subcritical  evaporation  behavior. 


(3) 


where  D  denotes  the  mean  self  diffusion  of  the  system  and  is  the  velocity  of  the  i*^  particle. 
This  mean  value  is  obtained  by  summing  the  velocity  autocorrelation  function  of  each  of  the 
particles  over  the  total  number  N. 

The  shear  viscosity  coefficient  measures  the  resistance  of  the  fluid  to  a  shearing  force.  For 
the  shear  viscosity  the  integration  is  over  the  non-diagonal  terms  of  the  stress  tensor  and  it  is  a 
collective  property  of  the  fluid.  The  Green-Kubo  formula  is  given  by, 

i  ('« 

0 

where  an  off-diagonal  term  of  the  stress  tensor  is  given  by 

t=l  ^  t=l  ;=l 

Here  T  and  V  denote  the  temperature  and  volume  of  the  system,  respectively,  and  m  denotes  the 
mass  of  each  particle,  v^,  rij  and  F(rij)  denote  the  x-direction  velocity,  position  vector  and  the 
potential,  respectively,  between  particles  i  and  j.  Statistical  precision  is  improved  by  averaging 
over  all  six  terms  that  result  from  the  stress  tensor: 


The  thermal  conductivity  coefficient  measures  the  transport  of  heat  in  a  system.  The 
correlation  function  is  obtained  from  the  heat  current  and  it  is  also  a  collective  property  of  the  entire 
fluid.  The  Green-Kubo  formula  is  given  by: 

^  0 

where  T*  denotes  an  arbitrary  component  of  the  heat  current  and  is  given  by 


where  i  denotes  a  unit  tensor.  Again,  statistical  precision  is  improved  by  averaging  over  all  three 
components: 

+  +  (9) 

These  functions  are  computed  during  the  simulation.  The  velocity  autocorrelation  is 
computed  outside  the  force  subroutine,  while  the  other  components  of  the  other  autocorrelation 
functions  are  calculated  within  it.  Once  the  correlation  functions  have  been  obtained,  standard 
numerical  integration  techniques  are  used  to  obtain  the  transport  coefficients.  The  statistical 
precision  of  the  self-diffusion  coefficient  is  improved  by  averaging  over  all  the  particles  in  the 
system.  The  viscosity  and  thermal  conductivity  values  are  improved  by  averaging  over  all  the  six 
terms  from  the  stress  tensor,  and  the  three  terms  of  the  heat  current,  respectively. 

The  MD  program  uses  the  effective  Lennard-Jones  potential,  with  system  sizes  of  256 
molecules,  and  simulations  of  500,000  timesteps  for  the  transport  coefficients  computation  and 
50,000  timesteps  for  pressures.  A  typical  timestep  size  of  2  fs  is  used  for  most  of  the  simulations. 
The  code  also  uses  linked  cell  lists  for  efficient  sorting  of  molecules,  periodic  boundary  conditions 
and  a  modified  velocity  Verlet  algorithm  for  particle  displacement. 


Simulations  have  been  carried  out  on  pure  argon,  nitrogen,  and  oxygen  at  various 
supercritical  conditions,  with  self-diffusion  coefficients,  shear  viscosity  coefficients,  thermal 
conductivity  coefficients,  and  pressures  computed  for  most  of  the  conditions.  Results  have  been 
compared  to  the  National  Institute  of  Standards  and  Technology  (NIST)  [9]  values.  Figure  6 
shows  diffusion  coefficients  for  oxygen.  It  can  be  observed  that  the  coefficient  decreases  with 
temperature  at  3  MPa,  but  increases  at  supercritical  pressures.  Figures  7  and  8  show  shear 
viscosity  and  thermal  conductivity.  The  MD  results  compare  well  with  the  NIST  values,  especially 
at  the  higher  temperatures.  The  poorer  match  at  160  K  for  the  10  and  15  MPa  cases  may  be  due  to 
the  use  of  an  insufficient  sample  size  as  the  critical  temperature  is  approached.  Figure  9  shows 
calculated  pressures  as  a  function  of  temperature  for  various  densities.  For  most  of  the  cases  there 
is  excellent  agreement  between  the  MD  and  NIST  values. 

Three-Bodv  Forces.  MD  simulations  normally  use  two-body  (pair-wise  additive) 
potentials  such  as  the  Lennard-Jones  12-6  potential  to  calculate  interatomic  forces.  This  is  done 
because  good  agreement  with  experimental  results  can  be  achieved  with  two-body  potentials  and 
because  calculations  using  three-body  potentials  become  extremely  computation^ly  intensive. 
However  concern  has  been  expressed  that  the  high  densities  associated  with  supercritical 
conditions  could  lead  to  errors  when  using  only  two-body  potentials.  Thus  an  argon  evaporation 
model  using  the  classical  three-body  Axilrod-Tellerpotential  [10]  has  been  constructed  in  order  to 
examine  any  difference  in  evaporation  behavior  under  supercritical  conditions  between  two-  and 
three-body  potentials  [11].  Little  difference  was  detected  between  the  two-body  and  three-body 
simulations. 


Fig.  6  Oxygen  seif-diffusion  coefficients  as  a  Fig.  7  Oxygen  shear  viscosity  as  a  function  of 

function  of  temperature  for  various  temperature  for  various  pressures 

pressures  calculated  via  MD  comparing  MD  and  NIST  calculated  values 


Fig.  8  Oxygen  thermal  conductivity  as  a  function 
of  temperature  for  various  pressures 
comparing  MD  and  NIST  calculated  values 


Fig.  9  Oxygen  pressure  as  a  function  of 
temperature  for  various  densities 
comparing  MD  and  NIST  calculated 
values 
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